rm(list=ls(all=TRUE))

library(ggplot2)
library(grid)
library(gridExtra)
library(splines)

setwd("YOUR-DIRECTORY")
dir.create("Figures")

##############
## FIGURE 1 ##
##############

dat1 <- read.csv(file="z3.PerMission")

dat1 <- dat1[order(dat1$abuses),]
dat1$mission <- as.character(dat1$mission)
dat1$mission <- factor(dat1$mission, as.character(dat1$mission))

gg1 <- ggplot(dat1, aes(x=mission, y=abuses)) +
    geom_bar(stat="identity", width=.65, color="#2b8cbe", fill="#2b8cbe") +
    ylab("Number of abuses") +
    xlab("UN mission title") +
    theme(
        axis.text.y=element_text(size=10, color="black"),
        axis.text.x=element_text(size=9, angle=45, hjust=1),
        axis.title.x = element_text(size=14, margin = margin(t=35, r=0, b=0, l=0)),
        axis.title.y = element_text(size=14),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm")
    )

dat2 <- read.csv(file="z4.PerCountry")
dat2$names <- c("Western Sahara", "Haiti", "Guatemala", "Yugoslavia",
                "Bosnia-Herzegovina", "Cyprus", "Georgia", "Ivory Coast",
                "Liberia", "Sierra Leone", "Central African Repub.", "Chad",
                "Dem. Repub. of Congo", "Burundi", "Eritrea", "Angola",
                "Sudan", "Iraq", "Syria", "Lebanon", "Israel", "East Timor")

dat2 <- dat2[order(dat2$abuses),]

dat2$names <- factor(dat2$names, as.character(dat2$names))

gg2 <- ggplot(dat2, aes(x=names, y=abuses)) +
    geom_bar(stat="identity", width=.65, color="#2b8cbe", fill="#2b8cbe") +
    ylab("Number of abuses") +
    xlab("Host country") +
    theme(
        axis.text.y=element_text(size=10, color="black"),
        axis.text.x=element_text(size=9, angle=45, hjust=1),
        axis.title.x = element_text(size=14),
        axis.title.y = element_text(size=14),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm")
    )

pdf("Figures/Fig1a.pdf", width=7, height=5, paper="special")
par(mar=c(0,0,0,0))
gg1
dev.off()

pdf("Figures/Fig1b.pdf", width=7, height=5, paper="special")
par(mar=c(0,0,0,0))
gg2
dev.off()

rm(list=ls(all=TRUE))


#########################
## FIGURES 2, 3, AND 4 ##
#########################

## Make a list of labels, one element for each graph
labs <- list(
    c("Contributor rule of law", "Contributor press freedom",
      "Host rule of law", "Host press freedom"),
    c("Contributor maternal mortality","Contributor economic inequality",
      "Contributor secondary enrollment"),
    c("Contributor IHL legislation","Contributor IHL enforcement",
      "Contributor conscription")
)

## Use standard grids, not facets, so each graph has its own axes
for ( i in 1:3 ) {
    m <- read.csv(paste("FromStata/margins_model", i, ".csv", sep=""))
    x <- 1:nrow(m)
    for (j in seq(1,ncol(m),2)) {
        if (j == 1 ) iter <- 1
        ## Use splines to get a smooth set of y points
        y <- predict(interpSpline(x,m[,j]))$y
        ## Do the same for x points and CIs
        xx <- predict(interpSpline(x,x))$x
        uci <- predict(interpSpline(x,m[,j]+(1.64*m[,j+1])))$y
        lci <- predict(interpSpline(x,m[,j]-(1.64*m[,j+1])))$y
        ## Turn it all into a data frame for ggplot
        use <- data.frame(
            mn = x,
            ymn = m[,j],
            yhi = m[,j] + ( 1.64*m[,j+1] ),
            ylo = m[,j] - ( 1.64*m[,j+1] )
        )
        gg.obj <- ggplot(use, aes(x=mn,y=ymn)) +
            geom_line(colour="#2c7bb6") +
            geom_point(colour="#2c7bb6", shape=1, size=1) +
            geom_ribbon(aes(ymin=ylo,ymax=yhi), alpha=0.2) +
            ylab("Predicted annual number of abuses") +
            xlab("Percentile") +
            scale_x_continuous(
                labels=paste(seq(10,90,10),"th",sep=""),
                breaks=use$mn
            ) +
            ggtitle(labs[[i]][iter]) +
            theme(
                axis.title.y = element_text(vjust=1, size=9),
                axis.title.x = element_text(vjust=-0.25),
                plot.title = element_text(size = 10)
            )
        assign(paste("gg", iter, sep=""),  ggplot_gtable(ggplot_build(gg.obj)))
        rm(gg.obj,y,xx,uci,lci,use)
        iter <- iter+1
    }
    if ( i == 1 ) {
    	plots <- list(gg1,gg2,gg3,gg4)
    	pdf("Figures/Fig2.pdf", width=10, height=6, paper="special")
    	do.call(grid.arrange, c(plots, list(ncol=2, widths=c(1,1))))
        dev.off()
        rm(gg1,gg2,gg3,gg4)
    }
    if ( i == 2 ) {
    	plots <- list(gg1,gg2,gg3)
    	pdf("Figures/Fig3.pdf", width=10, height=3.5, paper="special")
    	do.call(grid.arrange, c(plots, list(ncol=3, widths=c(1,1,1))))
        dev.off()
        rm(gg1,gg2,gg3)
    }
    if ( i == 3 ) {
    	plots <- list(gg1,gg2,gg3)
    	pdf("Figures/Fig4.pdf", width=10, height=3.5, paper="special")
    	do.call(grid.arrange, c(plots, list(ncol=3, widths=c(1,1,1))))
        dev.off()
        rm(gg1,gg2,gg3)
    }
}


##############
## FIGURE 5 ##
##############

rm(list=ls(all=TRUE))

dat <- read.csv("FromStata/Forecasts.csv")

## Grab the appropriate forecasts
dat1 <- rbind(
    dat[dat$year==2014 & dat$mission=="ForecastUKR",],
    dat[dat$year==2012 & dat$mission=="ForecastSYR",],
    dat[dat$year==2006 & dat$mission=="ForecastSPN",],
    dat[dat$year==2013 & dat$mission=="ForecastTUR",],
    dat[dat$year==2015 & dat$mission=="ForecastCOL",]
)

## Put everything into a data frame
use <- data.frame(
    predict = c(dat1$PredictBestUNOMIG,
                dat1$PredictAvg,
                dat1$PredictWorstCONGO),
    mission = rep(dat1$mission,3),
    scenario = c(rep("Best",5),rep("Avg",5),rep("Worst",5)),
    x = rep(c(2,5,1,4,3),3)
)
use$scenario <- factor(use$scenario, c("Best","Avg","Worst"))

gg.out <- ggplot(use, aes(x=x, y=predict, fill=scenario)) +
    geom_bar(stat="identity", position="dodge", width=.65) +
    scale_fill_manual("Scenario", values=c("#feb24c","#fc4e2a","#bd0026"),
                      labels=c("Best case", "Average case", "Worst case")) +
    scale_x_continuous(breaks=c(1,2,3,4,5),
                       labels=c("Spain\n2006","Ukraine\n2014","Colombia\n2015",
                                "Turkey\n2013","Syria\n2012")) +
    xlab("Hypothetical mission host country") +
    ylab("Predicted number of abuses") +
    theme(
        axis.text.x=element_text(size=9),
        legend.position=c(.2,.8)
    )

pdf("Figures/Fig5.pdf", width=5, height=4, paper="special")
par(mar=c(0,0,0,0))
gg.out
dev.off()

rm(dat1,gg.out,use)

##############
## FIGURE 6 ##
##############

years <- 2002:2015

dat2 <- rbind(
    dat[dat$mission=="ForecastTUR" & dat$year >= years[1] & dat$year <= years[length(years)],]
)
dat2 <- dat2[,c("year","PredictWorstCONGO")]

gg.out1 <- ggplot(dat2, aes(x=year, y=PredictWorstCONGO)) +
    geom_bar(stat="identity", width =.65, fill="#bd0026") +
    xlab("Year") +
    ylab("Predicted abuses (worst case)") +
    scale_x_continuous(labels=dat2$year, breaks=dat2$year) +
    theme(
        axis.text.x = element_text(angle = 45, hjust = 1)
    ) +
    ggtitle("Out-of-sample predictions")

dat3 <- data.frame(
    vals = c((dat[dat$mission=="ForecastTUR" & dat$year >= years[1] & dat$year <= years[length(years)], "pressFreedomHC"] - dat[dat$year==2002 & dat$mission=="ForecastTUR",]$pressFreedomHC) / 
    dat[dat$year==2002 & dat$mission=="ForecastTUR",]$pressFreedomHC * 100,
    (dat[dat$mission=="ForecastTUR" & dat$year >= years[1] & dat$year <= years[length(years)], "ruleoflawHC"] - dat[dat$year==2002 & dat$mission=="ForecastTUR",]$ruleoflawHC) /
    dat[dat$year==2002 & dat$mission=="ForecastTUR",]$ruleoflawHC * 100
    ),
    vars = c(
        rep("Press freedoms",nrow(dat2)),
        rep("Rule of law",nrow(dat2))
    ),
    yr = rep(dat2$year,2)
)

gg.out2 <- ggplot(dat3, aes(x=yr, y=vals, color=vars)) +
    geom_line() +
    geom_point(aes(shape=vars)) +
    ylim(min(dat3$vals),max(dat3$vals)) +
    xlab("Year") +
    ylab("Percentage change since 2002") +
    scale_color_manual(values=c("#3288bd", "#99d594")) +
    scale_x_continuous(labels=dat2$year, breaks=dat2$year) +
    theme(
    	legend.position=c(0.2,0.25),
    	legend.title=element_blank(),
    	axis.text.x = element_text(angle = 45, hjust = 1)
    ) +
    ggtitle("Observed attributes")

plots <- list(gg.out1,gg.out2)
pdf("Figures/Fig6.pdf", width=10, height=3.5, paper="special")
par(mar=c(0,0,0,0))
do.call(grid.arrange, c(plots, list(ncol=2, widths=c(1,1))))
dev.off()

rm(dat2,dat3,gg.out1,gg.out2,plots,years)


##############
## FIGURE 7 ##
##############

years <- 2002:2015

dat2 <- rbind(
    dat[dat$mission=="ForecastCOL" & dat$year >= years[1] & dat$year <= years[length(years)],]
)
dat2 <- dat2[,c("year","PredictWorstCONGO")]

gg.out1 <- ggplot(dat2, aes(x=year, y=PredictWorstCONGO)) +
    geom_bar(stat="identity", width =.65, fill="#bd0026") +
    xlab("Year") +
    ylab("Predicted abuses (worst case)") +
    scale_x_continuous(labels=dat2$year, breaks=dat2$year) +
    theme(
        axis.text.x = element_text(angle = 45, hjust = 1)
    ) +
    ggtitle("Out-of-sample predictions")

dat3 <- data.frame(
    vals = c(
    (dat[dat$mission=="ForecastCOL" & dat$year >= years[1] & dat$year <= years[length(years)], "pressFreedomHC"] - dat[dat$year==2002 & dat$mission=="ForecastCOL",]$pressFreedomHC) / 
    dat[dat$year==2002 & dat$mission=="ForecastCOL",]$pressFreedomHC * 100,
    (dat[dat$mission=="ForecastCOL" & dat$year >= years[1] & dat$year <= years[length(years)], "ruleoflawHC"] - dat[dat$year==2002 & dat$mission=="ForecastCOL",]$ruleoflawHC) /
    dat[dat$year==2002 & dat$mission=="ForecastCOL",]$ruleoflawHC * 100
    ),
    vars = c(
        rep("Press freedoms",nrow(dat2)),
        rep("Rule of law",nrow(dat2))
    ),
    yr = rep(dat2$year,2)
)

gg.out2 <- ggplot(dat3, aes(x=yr, y=vals, color=vars)) +
    geom_line() +
    geom_point(aes(shape=vars)) +
    ylim(min(dat3$vals),max(dat3$vals)) +
    xlab("Year") +
    ylab("Percentage change since 2002") +
    scale_color_manual(values=c("#3288bd", "#99d594")) +
    scale_x_continuous(labels=dat2$year, breaks=dat2$year) +
    theme(
    	legend.position=c(0.2,0.85),
    	legend.title=element_blank(),
    	axis.text.x = element_text(angle = 45, hjust = 1)
    ) +
    ggtitle("Observed attributes")

plots <- list(gg.out1,gg.out2)
pdf("Figures/Fig7.pdf", width=10, height=3.5, paper="special")
par(mar=c(0,0,0,0))
do.call(grid.arrange, c(plots, list(ncol=2, widths=c(1,1))))
dev.off()

rm(list=ls(all=TRUE))


